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ABSTRACT 

The  one-dimensional  roughness  of  an  under-ice  profile  of  elliptical 
bosses  is  modeled  in  the  time  domain  by  a  shot  noise  process  of  elliptical 
pulses  of  random  amplitude,  duration,  and  time  of  occurrence.  A  sample 
realization  of  8000  data  points  is  generated  and  plotted  for  visual  comparison 
with  experimental  under-ice  data.  Also,  theoretical  and  simulation  results 
for  the  power  density  spectrum,  the  autocorrelation  function,  the 
characteristic  function,  the  cumulative  distribution  function,  and  the 
probability  density  function  of  the  shot  noise  process  are  plotted  and 
compared. 
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INTRODUCTION 

The  under-ice  profile  has  been  observed  to  appear  like  a  random 
collection  of  superposed  elliptical  bosses,  each  of  random  amplitude,  length, 
and  location.  An  analogous  model  in  the  time  domain  is  shot  noise  composed  of 
overlapping  pulses  of  random  amplitude,  duration,  and  time  of  occurrence. 
Accordingly,  we  have  generated  a  sample  realization  of  a  shot  noise  process 
for  visual  comparison  with  experimental  under-ice  data,  and  for  possible 
corroboration  of  this  model.  The  particular  realization  generated  has  8000 
data  points,  although  the  number  of  effectively-independent  samples  is  far 
fewer,  as  will  be  demonstrated. 

A  number  of  analytical  results  for  shot  noise  have  been  derived  in  the 
past  [1];  however,  they  did  not  cover  the  case  of  random  duration  modulation. 
We  have  extended  the  analyses  to  include  random  durations  (as  well  as  random 
amplitudes  and  random  time  occurrences)  and  evaluated  the  spectrum  of  the  shot 
noise  process,  as  well  as  the  autocorrelation  function  and  the  first-order 
characteristic  function  of  the  instantaneous  amplitude.  From  the  latter,  the 
first-order  probability  density  function  and  cumulative  distribution  function 
of  shot  noise  have  been  evaluated  via  a  generalized  Laguerre  expansion 
employing  32  cumulants  or  moments.  Comparisons  of  all  these  theoretical 
results  with  the  corresponding  sample  quantities,  obtained  from  the  8000  data 
point  realization  above,  reveal  excellent  agreement. 
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A  REALIZATION  OF  A  SHOT-NOISE  PROCESS 

Shot  noise  is  characterized  by  a  superposition  of  pulses,  each  located 
independently  and  uniformly  on  the  time  scale.  A  sample  pulse  is  illustrated 
in  figure  1.  The  time  of  occurrence  t^  (center  of  symmetrical  pulse,  for 


example)  is  uniformly  distributed  in  time  t,  with  an  average  number  of  pulses 
per  second,  v.  The  amplitude  ak  and  half-duration 
pulse  are  also  all  independent  and  are  each  identically  randomly-distributed 
with  arbitrary  probability  density  functions.  Finally  the  fundamental  pulse 
shape  F  in  figure  1  is  arbitrary. 


of  an  individual 


A  realization  of  shot  noise  is  given  by 


I(t)  = 


(1) 


where  the  summation  extends  over  all  k.  The  particular  data  we  generate  here 
employs  the  following  example;  unsealed  pulse  shape  F  is  circular: 


F(x) 


(2) 
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This  pulse  is  continuous;  however,  it  has  cusps  (infinite  slope)  at  x  =  *1. 
The  reason  for  this  selection  will  become  apparent  when  we  discuss  the 
spectrum  of  shot  noise  process  (1). 


The  amplitude  probability  density  function  for  random  variable  a^  is 
Rayleigh, 

'  .21 


P (a )  =  “f  exp  )U(a)  , 


(3) 


and  the  duration  probability  density  function  for  random  variable^  is  also 
Rayleigh, 


PU)  =  4  exP  (4?  )  * 

%  \2V 


(4) 


Here,  step  function 


U(x)  * 


1  for  x  >  0 


(5) 


0  for  x  <  0  J 

The  mean  values  of  random  variables  a^  and  are  given  respectively  by 


a  "  ak  =  (f)CTa’  J  =  h  =  ff 


V 


(6) 


in  terms  of  the  parameters  afl  and  a  of  probability  density  functions 
(3)  and  (4).  Alternatively,  the  mean  square  values  are  given  by 


T  T  o  2  77  IT 


a  =  a 


-  2oa  *  A' =4  -  2o; 


* 


(7) 


Three  typical  component  pulses  are  depicted  in  figure  2,  and  can  range 


from  circular  through  various  elongated  elliptical  shapes.  The  total  length 
of  an  individual  pulse  is  =  2^.  An  important  parameter  of  this  time- 
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^3 


limited  pulse  shape  in  figures  1  and  2  is  the  (dimensionless)  overlap  factor 


(8) 


This  is  the  average  number  of  pulses  that  are  overlapping  at  any  one  instant 
of  time,  and  is  a  partial  measure  of  the  applicability  of  the  central  limit 
theorem.  A  more  meaningful  measure  are  the  cumulants;  for  probability  density 
function  (3)  and  pulse  shape  (2),  the  normalized  third  and  fourth  cumulants  are 


1.017 


and 


1.2 

A* 


(9) 


respectively.  In  the  sample  realization  generated  here,  the  overlap  factor  in 
(8)  was  6.2,  leading  to  normalized  cumulant  values  in  (9)  of  .58  and  .39, 
respectively.  Since  a  Gaussian  probability  density  function  would  lead  to 
zero  cumulants  above  second-order,  the  shot  noise  realization  dealt  with  here 
is  distinctly  non-Gaussian. 
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In  the  three  parts  af  figure  3,  a  realization  of  shot  noise  model  (1)  is 
given  for  parameter  values 

0a  =  1  sec,  =20  sec,  v  =  .124  pulses/sec.  (10) 

The  waveform  (1)  is  sampled  at  unit  time  increments  and  connected  by  straight 

lines;  thus  the  initial  100  data  points  illustrated  in  figure  3A  have  a  jagged 
appearance  for  those  component  pulses  with  small  as  for  example  at  time 
instants  67-68.  The  larger  duration  pulses,  like  the  one  centered  at  t  =  29, 
have  a  smoother  appearance. 

In  figure  3B,  the  initial  1000  data  points  illustrate  the  very  erratic 
character  of  shot  noise;  the  waveform  consists  of  some  very  sharp  spiky  pulses 
and  other  broader  smooth  components.  The  appearance  of  a  downward  trend  in 
these  1000  data  points  is  erased  when  the  entire  8000  data  point  sequence  is 
viewed  in  figure  3C.  The  possibility  of  shot  noise  process  (1)  reaching  a 

zero  value  (when  no  pulses  overlap)  is  confirmed  by  the  waveform  values  near 

t  =  6400. 


/ 
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CORRELATION  AND  SPECTRUM  OF  SHOT  NOISE  PROCESS 


The  derivations  of  the  correlation  and  spectrum  of  the  shot  noise  process 
(1)  are  given  in  appendix  A;  from  (A-12),  we  have,  in  general,  the  correlation 
function  at  delay  T, 


Rj(T)  =  v  7  f  d#  p  4>(TU)  +  Idc  , 
where  the  dc  component  of  I(t)  is,  from  (A-13), 

!dc  =vaj  Jdx  F  (x )  , 


(ID 


(12) 


and 


<(>(y)  =  J  dx  F(x)  F(x-y) 

is  the  (aperiodic)  correlation  of  an  individual  pulse  F. 
over  the  range  of  non-zero  integrand.) 


(13) 

(All  integrals  are 


Also,  from  (A-16),  the  general  spectrum  of  process  I(t)  is,  at  frequency 

f. 


Gj(f)  =  v  7  J  pte)  A2  ]stff)|2  +  Idc  S(f)  »  (14) 


where 


S(f)  =  Jdx  exp(-i2irfx)  F(x)  (15) 

is  the  voltage  density  spectrum  (Fourier  transform)  of  pulse  F.  Thus 
|s(f)p  is  the  energy  density  spectrum  corresponding  to  pulse  F. 
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It  should  be  observed  that  the  entire  probability  density  function  pf£) 
of  half-duration  random  variable  is  required  in  order  to  evaluate  the 

correlation  or  spectrum  of  shot  noise.  However,  only  the  first  two  moments,  "a 

— n 

and  a  ,  are  required  known  about  probability  density  function  p (a )  of 
amplitude  random  variable  a^.  The  only  way  that  the  dc  term  Idc  can  be 
zero  is  if  random  variable  a^  has  zero  mean  fa  =  0),  or  if  pulse  F  has  zero 
area  (S(0)  =  0). 

Example 


The  example  of  interest  here  was  given  earlier  in  (2)  and  (4),  namely  a 
circular  pulse  F  and  a  Rayleigh  probability  density  function  for  random 
variable^.  The  spectrum  6j(f)  in  (14)  is  evaluated  in  (A-17)  through 
(A-22),  with  the  results 


Jl(2irf> 

S(f)  =  ’  S(0)  -  £ 


dc 


/.f  - 
-{? j  va°- 


i.  • 


r  m  „2  7  22exp(-z)  1,(2)  ,3  2  2  2 

6j  (f )  =  2ir  v  a  - - -  {j)  V  3  ^  ®  f 

with  z  =  (2tt<^  f  )2.  (16) 

The  asymptotic  behavior  of  spectrum  (16)  is  [2,  eq.  9.7.1] 

~3  , 

Gt(f)^  f  as  f -*+•«.  (17) 


(2ir) 


5 


14 


TM  No.  841208 


That  is,  the  spectrum  decays  at  a  -30  dB/decade  rate  at  large  frequencies; 
this  is  due  to  the  square  root  singularities  at  x  =  ±1  of  pulse  F  given  in 
(2).  This  decay  rate  has  been  observed  in  some  spectral  analyses  of  under -ice 
profiles,  and  was  one  of  the  reasons  for  choosing  the  specific  circular  pulse 
in  (2)  for  this  investigation. 


The  spectrum  in  (16)  is  plotted  in  figure  4,  for  the  choice  of  parameters 
earlier  in  (10),  as  a  dashed  line,  normalized  to  0  dB  at  f  =  0.  Superposed  is 
a  1  inear -predictive  spectral  analysis  result  with  predictive  order  10,  for  the 
8000  data  points  of  figure  3C.  The  two  results  are  in  excellent  agreement, 
even  at  the  -50  dB  level,  with  the  inevitable  3  dB  aliasing  effect  at  the 
Nyquist  frequency,  as  indicated. 


The  correlation  Rj(T)  in  (11)  is  evaluated  in  (A-23)  through  (A-33), 
for  the  example  (2)  and  (4),  with  the  result 


Rj(r)  =  |{2irfl  s  exp(-s)  [  (l+4s )  Kx(s)  -  (3+4s)  KQ(s )]  + 


+ 


(18) 


This  quantity,  exclusive  of  the  I^c  term,  and  normalized  at  the  origin,  is 
plotted  in  figure  5  as  a  dashed  line,  for  delays  (lags)"C  up  to  100.  It  is 
seen  to  decay  monotonical ly  to  zero  as  T  increases,  and  reach  its  1/e  value  at 
approximately  T=  30. 
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The  remaining  solid  curve  on  figure  5  is  the  normalized  sample 
autocorrelation  function  of  the  8000  data  point  sequence  in  figure  3C,  where 
the  sample  mean  was  subtracted  from  the  given  data.  The  agreement  with 
theoretical  result  (18)  is  excellent.  The  dotted  horizontal  lines  at  ±2o  in 
figure  5  are  the  ±2  sigma  values  of  the  correlation  estimate  at  delays  where 
the  true  correlation  is  presumed  zero;  the  details  of  this  analysis  are  given 
in  appendix  B. 


This  procedure  is  duplicated  in  figure  6,  where  the  correlation  function 
estimate  out  to  lag  TT=  1000  is  plotted.  The  drifting  of  the  estimate  outside 
the  ±2ct  limits  (at  T=  470  and  820)  is  consistent  with  an  occasional  excursion 
of  a  random  variable  outside  its  ±2o  range.  The  correlation  estimate  (used 
for  figures  5  and  6)  at  time  separation  k  is 


N 

Rl  -  i  .SE.  x  x  .  for  k  >  0  , 
K  n  n=k+l  n  n-K 


(19) 


where  M,  is  the  available  data  in  figure  3C,  with  its  sample  mean  removed. 
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AMPLITUDE  STATISTICS  OF  SHOT  NOISE 


The  first-order  characteristic  function  of  shot  noise  process  I(t)  is 
derived  in  appendix  C;  it  is  given  by  (C-9)  as 

fj(5)  =  exp  [v  jJdx{fa5F(x)]-l}]  .  (20) 

Here  f.  is  the  first-order  characteristic  function  of  amplitude  random 

a 

variable  a^.  Observe  that  the  probability  density  function  p(^)  of 
duration  ^  is  irrelevant  to  characteristic  function  f  j,  except  for  its 
mean  J;  this  is  in  contrast  to  the  spectrum  and  correlation  results  in  (11) 
and  (14),  where  p(a)  was  irrelevant  except  for  parameters  a  and  a2.  (For 
/^  =  1  for  all  k,  (20)  reduces  to  a  simplified  version  of  [1,  eq.  1.5-4].) 


The  characteristic  function  of  the  amplitude  random  variable  a^  can  be 
expanded  in  terms  of  its  moments 


ua(n)  a  an  =  J" da  an  p(a)  for  n  2  0  » 


(21) 


according  to 


f  (?)  -  ^L».(n)d ?>"/»; 


a  5  fT'a 
n=sU 


(22) 


This  result  is  useful  if  the  ^n  of  (20)  is  expanded  in  a  series  inj;  namely 


f  At)  =  v  ua(n)(i|)n 

n=l 


dx  Fn(x)/ni  , 


(23) 


giving  immediately  the  cumulants  of  I(t)  as 


/^(n)  =  v  X  wa(n)  J  dx  F°(x)  for  n  2  1  • 


(24) 
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That  is,  the  n-th  cumulant  of  I(t)  is  proportional  to  the  n-th  moment  of 
random  variable  a^  as  well  as  the  n-th  “moment"  of  pulse  F.  (For  ^  =  1 
for  all  k,  (24)  reduces  to  [1,  eq.  1.5-2].) 


The  normalized  cumulant  of  I(t)  is 

Yr(n)  - 


%(n)  •  1  ua(n)  j 

dx  Fn(x) 

"A  *  /  _£-»  ' 
CCj(2)]  (vjO 

[“a<2)  . 

f  dx  F2(x)J 

!*/* 

(25) 


In  particular,  the  coefficients  of  skewness  and  excess  [3,  pp.  184  and  187] 
are 


rI 


(3)  = 


(vi) 


Vi 


va(3)  , 

dx  F3(x) 

- 

r  2  l3/z 

ua(2) 

a  J 

L. 

dx  F^(x) 

(26) 


and 


i  ua(4 )  f dx  F4(x ) 

Yj  (4 )  -jl - - zr 

v  l  ua(2)  j  dx  F2(x)J 


(27) 


These  quantities  are  very  important  measures  of  the  approach  of  I(t)  to  a 
Gaussian  process;  if  vJL  is  very  large,  the  normalized  cumulants  Yj(n)  are 
all  substantially  zero  for  n  >_  3,  meaning  that  I(t)  is  nearly  Gaussian.  Thus 
although  probability  density  function  p(£)  is  not  directly  relevant  to  the 
probability  density  function  or  characteristic  function  (20)  of  I(t),  the 
exact  probability  density  function  of  I(t)  is  critically  dependent  on  the 
mean  Ji  through  the  dimensionless  parameter  vJL.  More  precisely,  (26)  and  (27) 
are  the  critical  quantities;  see  also  [1,  eq.  1.6-3]. 
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If  either  the  third  moment  of  random  variable  a,,  is  zero,  or  if  the 
third  moment  of  pulse  F  is  zero,  then  yj(3)  =  0.  In  that  case,  yj(4)  is 
the  most  important  statistic  measuring  the  applicability  of  the  central  limit 
theorem;  yj(4)  can  never  be  zero  for  shot  noise,  since  neither  the  fourth 
moment  of  random  variable  a^  or  pulse  F  can  be  zero  (except  in  a  trivial 
case). 


The  first  moment  of  shot  noise  I(t)  is  the  mean 

F(X)  (28) 

and  has  already  been  encountered  in  (12).  It  can  be  zero  only  if  the  first 
moment  of  random  variable  ak  or  of  pulse  F  is  zero. 

Example 


Numerous  cases  have  been  considered  in  appendix  C;  in  the  main  body  here, 
we  limit  attention  to  example  (2)  and  (3)  presented  earlier.  We  find 


M 

»a(n)  -  Z7r(j  +  l)  for  n  >  0  , 


r  2n+1  r2  (7  *  A 

jdx  F"(x) - 


for  n  >  0  . 


(29) 


Then  (26)  and  (27)  yield  result  (9)  quoted  earlier. 


The  realization  of  shot  noise  process  I (t )  in  figure  3C  employed  the 
parameters  in  (10).  The  sample  cumulative  distribution  function  of  these  8000 
data  points  is  depicted  in  figure  7,  on  a  normal  probability  ordinate;  thus  a 
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truly  Gaussian  random  variable  would  have  the  straight  line  character 
indicated.  The  significant  deviation  of  the’sample  cumulative  distribution 
function  from  the  Gaussian  line  is  due  to  the  small  value  of  the  overlap 
factor  in  (8),  namely 


ITk  v  =  2  Jk  v  =  6.2  .  (30) 

The  moments  in  (29)  are  all  positive  and  are  easily  numerically  evaluated 
via  recursion;  hence  the  cumulants  in  (24)  can  be  accurately  evaluated  for 
high-order  n.  When  these  cumulants  are  employed  in  a  generalized  Laguerre 
expansion  of  the  cumulative  distribution  function  of  I(t),  using  32  moments  of 
(29),  the  solid  curve  in  figure  8  is  obtained.  The  sample  cumulative 
distribution  function  of  figure  7  is  duplicated  here,  although  the  abscissa  is 
scaled  differently.  The  agreement  between  theory  and  experiment  in  figure  8 
is  excellent,  considering  the  fact  that  we  only  have  about  8000/30  =  270 
effectively  independent  samples  of  I(t)  in  figure  3C;  the  denominator  factor 
of  30  here  is  the  effective  correlation  duration,  previously  identified  in 
figure  5  at  the  1/e  point. 

Finally,  when  the  same  32  moments  are  used  in  a  generalized  Laguerre 
expansion  of  the  probability  density  function  of  I(t),  the  result  in  figure  9 
is  obtained.  The  small  bump  near  the  origin  is  real  and  accurate;  it  and  the 
non-symmetric  tails  of  the  probability  density  function  confirm  the  distinctly 
non-Gaussian  character  of  I(t).  The  method  for  the  determination  of  the 
cumulative  distribution  function  and  probability  density  function  in  figures  8 
and  9  will  be  presented  in  a  NUSC  Technical  Report  [4]  by  the  author;  the 
programs  are  listed  here  in  appendix  0,  along  with  an  example  of  the  sequence 
of  Laguerre  coefficients. 
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APPENDIX  A.  DERIVATION  OF  SPECTRUM  AND  CORRELATION 


The  method  employed  below  follows  that  given  by  Rice  [1,  sections  1.4  and 
1.5]  rather  closely.  We  generalize  [1,  eq.  1.3-1]  to  the  current  form 
introduced  in  (1): 


K  ■ 

IK(t)  =2  ak 
K  k=l  k 


(A-l ) 


•  M’ 

presumed  number  of  pulses  to  occur  in  a  large  time  interval  T,  and  ak  is  a 
random  amplitude  as  in  [1,  eq.  1.5-1];  but  random  duration  ^  is  new.  Then 
product 


where  [a^ 


{41  are  all  independent  random  variables.  K  is  the 


Holding  random  variables 


{a^  and  fixed  for 


now,  the  statistical 


average  of  ( A— 2 )  over  ^t^  is 


k^m 


K  K  K 

-  f  2  ak  4  ffr/'V  +  Zz  -S.  2.  ak  am  4  4n  s2(0)  »  (A_3) 

k=l  1  k=l  m=l 

kfm 
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where  T  is  an  arbitrary  large  (but  finite)  time  interval,  and 

<j>(y)  =  j*dx  F(x)  F(x-y)  (A-4) 

is  the  aperiodic  autocorrelation  of  pulse  F,  while 

S(f )  =  j'dx  exp(-i2irfx)  F(x)  (A-5) 

is  the  voltage  density  spectrum  of  F. 

The  remaining  averages  over  independent  random  variables  [a^  and 
in  (A  -3)  now  yield 


j  K  ?  [dje  pUU  tfciJL)  +  (K2-K)  [a  ji  S(0)]2,  (A-6 

where  p(j>)  is  the  probability  density  function  of  random  variable  /|<* 


) 


Now  K  is  itself  a  random  variable,  with  discrete  probability 
(in  an  interval  T)  of  [1,  eq.  1.1-3] 

-kilj-exp(-vT)  for  K  -  0,  1,  2 . 


(A— 7 ) 


There  then  follows  the  characteristic  function  of  random  variable  K  as 

fK(£)  -  exp(vT[exp(i£)  -  1]),  (A-8) 


with  series  expansion 


oc 


JU  fK(?)  =  VT  [exp(i?)-l]  =  vT  (if)"/n!  .  (A-9) 

n=l 


Thus  the  cumulants  of  random  variable  K  are  all  equal, 

%K(n)  =  vT  for  n  >  1  , 

giving  in  particular  the  first  two  moments 


( A— 1 0) 


K  =  vT,  K  =  vT (vT  +  1)  . 
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The  use  of  (A-ll),  to  perform  the  remaining  average  of  (A-6)  with  respect 
to  random  variable  K,  then  yields  the  correlation  function  of  the  shot  noise 
process  I(t): 

Rj(D  =  v  7  J"dj l  piX)^  ZIJ> )  +  [v  S(0)]2.  (A— 12 ) 

The  dc  component  of  I(t)  is 

Idc  =  v  *  X  S(°)  =  v  aj  J  dx  F(x)  .  (A-13 ) 

The  spectrum  of  I(t)  is  the  Fourier  transform  of  ( A— 12 ) : 

Gj(f )  -  v  7  j*  dJt  p(£)  X2fW  )  +  Ijc  S(f)  .  (A— 14 ) 

where 

$(f )  -  J  dr  exp(-i2irfT)  <J>(T)  = 

=  j  dr  exp(-i2irfT)  J  dx  F(x)  F(x-r)  =  js(f)J2,  (A-15) 
by  use  of  (A-4)  and  (A-5).  Thus  (A-14)  can  be  expressed  as 


( A— 16 ) 


(A— 17 ) 
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Then  from  (A-5)  and  [5,  eq.  3.752  2], 

dx  exp(-i2^fx}  (l-x  J= — -  • 

Substitution  of  (A-17)  and  (A-18)  in  the  integral  in  (A-16)  yields, 
[5,  eq.  6.633  2], 


where 


2  exp(-z)  Ix(z) 
z 


» 


z  =  (2TO^f)2. 

Then  the  spectrum  (A-16)  is  given  by 


?  -a-  o  2  exp(-z)  I.  (z)  ? 

Gj(f)  =  2ir  v  a2  o2  - - - - - +  I2C  g(f)  , 


where 


- 


*dc  =  V  3  *  S(0)  =  (f)  V  a 


by  means  of  (A-13),  (A-18),  and  (6). 

To  determine  the  correlation  of  shot  noise  process  I(t),  we 
first  the  continuous  portion  of  the  spectrum  in  (A-21): 


2  -£  2  2  exp(-4ir2ajj2f2)  1-^  (4ir2^2f 2 ) 

G  (f)  =  2tt  V  a  a  - yyv ; - 

'  4*¥f 


(A-18) 

y  use  of 


(A-19) 


(A-20) 


(A-21) 


(A-22 ) 


consider 


(A-23) 


30 


TM  No.  841208 


The  corresponding  correlation  is 


+0* 

RC(D  =  J~  df  exp(i2irfr)  Gc(f) 
-•© 


,2  a2  2 
=  4  IT  v  a  CT^ 


df  2  C0S(2irfT) 


exp(-4,r2op2f2)  1^  (4ir2o2f2 ) 

A  2  2.-2 
4ir  f 


=  2  ir  v  a  or 

/ 


/t  A  exp(-z)  ‘i(z> 

|dz  C0S(°PV  P5 


yfi  o 

=  2  H  V  a  a  exp(-s)  W  %  (2s) 

“y  >  » 


where  we  employed  (A-20)  and  [5,  eq.  6.755  2],  and  defined 

2 


s  = 


The  W-function  in  (A-24)  is  the  Whittaker  function  [2,  p.  505]. 
Now  by  [5,  eqs.  9.232  1  and  9.222  1],  we  have 


°0 

W  3  1(2s)  =  W  3  x(2s)  =  f  dt  exp(-2st)  t/a  (1+t) 

~2*  2  ~2'~2  o 


s/a 


[ 


=  -|exp(-s)  j  dt  exp(-2st) 


’/a  %  f/a 

t  (1+t)  t  (1+t) 


But  according  to  [5,  eq.  3.364' 3], 


f  dt  PXP{~2st^  =  exp(as)  KQ(as) 
“o  t  a(a+t)  2 


(A-24) 


( A— 25 ) 


(A-26) 


(A— 27 ) 
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Partial  differention  with  respect  to  a  then  yields 

06 

-  2*  J dt  6^P^~2S^  =  s  e>^P(as)  [KQ(as )  -  <1(as )] 


0  t  (a+t) 


and  (repeated) 


r  J  dtSf^-2  •■*(«) 


t  (a+t) 


K.  (as) 

2K0(as)  -  2Kx(as)  + 


(A— 28 ) 


(A-29 ) 


Here  we  used  [2,  eq.  9.6.28]  in  the  forms 

Mz) 

Kq(z)  =  -  K^(z),  Kj (z)  =  -  K0(z)  -  .  (A— 30 ) 

If  we  now  set  a  =  1  in  (A-28)  and  (A-29),  and  then  employ  these  results 
in  (A-26),  we  obtain 


W3  i(2s)  -  ^-fs  [(1+4s)K1(s)  -  (3+4s)K0(s)].  (A-31) 


Finally,  the  use  of  (A-31)  in  (A-24)  yields 


Rc(r)  =  |(2ir)/a  V  a2  ^  s  exp(-s)[(l+4s)K1(s)  -  (3+4s)KQ(s )].  (A-32) 

The  Fourier  transform  of  the  impulsive  part  of  the  spectrum  in  (A-21)  is  simply 
the  constant 


(A-33) 


which  must  be  added  to  RC(T)  in  (A-32)  to  obtain  Rj(r).  Here  s  is  given  by 


(A— 25 ) . 

As  "C-»0+,  there  follows  from  (A-32), 


lim  R(T)  =  |(2^)’/av  7  cl  .  ( A— 34 ) 

T->0¥  C  J  2 
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APPENDIX  B.  VARIANCE  OF  CORRELATION  ESTIMATE 


Let  the  available  data  be  £xn*^,  with  zero  mean  and  variance  a^: 


—  T  2 

x„  =  0,  x  =  a  for  1  <  n  <  N  . 
n  n  —  — 


(B-l) 


The  autocorrelation  estimate  at  delay  k  is  defined  here  as 


"k  -ff 


M  x  x  .  for  k  >  0 

N  n  n-k 


(B-2) 


At  delay  0,  the  mean  value  of  estimate  RQ  is 


n—  _  j-  ->  „2  2 

K„  =  tr  =  a  . 


T-  1 
^o  =  fT 


n=l 


( B— 3 ) 


We  now  want  to  evaluate  the  standard  deviation  of  estimate 
k  large  enough  that  xn  and  xfl_j<  are  statistically  independent, 
mean  value 


R^  at  delays 
We  have 


N 

^k  =  ~  ^  xn  xn-k  =  ^  ’ 
K  N  n=k+l  n  n  K 


( B— 4 ) 


using  the  independence  at  separation  k.  The  mean  square  value  of  estimate 
Rk  is 

'  N 

“k  -  7  VVWVT  •  ‘B-5> 

fr  m,n=k+l 

For  the  large  separation  values  k  of  interest  here,  the  only  statistical 
dependence  that  contributes  non-trivially  to  the  double  sum  is  the  following: 


33 


TM  No.  841208 


N 


N  m,n=k+l 


xm  xn  xm-k  xn-k 


4  9 

V  >  (m-n)  = 

N^  m,n=k+l 


4 


a 


5^  (N-k-JnJ  )  p2(n)? 
|n|  <N-k 


(N-k)J^  p2 (n ) 

n 


( B— 6 ) 


Here  p  is  the  correlation  coefficient  of  data  and  we  have  assumed  that 

N  is  moderately  larger  than  the  effective  correlation  length  of  p.  The  ratio 
of  the  standard  deviation  of  estimate  to  the  mean  value  at  k=0  is  then 
the  normalized  standard  deviation  at  separation  k: 


(B-7) 


Notice  that  no^  Gaussian  assumptions  on  data  ^x^  have  been  employed  in  this 
analysis;  however,  p(k)  is  essentially  zero  at  the  k  values  of  interest. 


As  an  example,  for  an  exponential  correlation  of  effective  length  K0, 
there  follows 

p^(n)  =  JE  exP  ^x  exP  =  Ke  »  ( B  —8 ) 

where  we  assume  that  Ke  is  moderately  laryer  than  unity.  Then  (B-7)  yields 

^C(N-k)Ke]^.  (B— 9  ) 

These  results  hold  only  for  those  values  of  k  where  p(k)  has  substantially 
gone  to  zero.  Larger  values  of  K0  lead  to  larger  relative  standard 
deviations;  th.s  is  consistent  with  the  fact  that  there  are  then  a  lesser 
number  of  effectively-independent  samples  in  the  limited  data  set  of  length  N. 
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For  the  8000  data  point  example  of  interest  here,  inspection  of  figure  5 
reveals  that  Kg  =  30.  Thus 

+  <11  +  (8000-1$  /n  irvl 

2ak  =  ~7T0  ’  (B-10) 

These  confidence  limits  are  superposed  as  dotted  lines  on  figures  5  and  6. 
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APPENOIX  C.  PROPERTIES  OF  CHARACTERISTIC  FUNCTION  OF  SHOT  NOISE 


Derivation  of  Characteristic  Function 


The  method  of  derivation  of  the  characteristic  function  of  I (t )  presented 
here  parallels  that  of  Rice  [1,  sections  1.4  and  1.5]  very  closely.  We 
generalize  [1,  eq.  1.3-1]  to 


(C-l) 


and  the  ensuing  discussion.  The  characteristic  function  of  an  individual 
component  in  (C-l)  is 


where 


h) 


are  all  independent  random  variables;  see  (A-l) 


r  /'t-tj.V) 

exp 

I1'*  Fhr)J 

(C-2 ) 


where  the  statistical  average  is  over  ak  ,  tk  ,  .  The  average  over 

tk  (for  fixed  ak,  is,  for  T  a  large  but  finite  time  interval 
[1.  P-  152], 


(C-3) 


37 


TM  No.  841208 


for  large  T,  where  we  have  used  the  fact  that 

F(x)-»  0  as  x-»  (C-4) 

Let  x  =  (t-Ijc)/JJ|<  in  (C-3)  to  get 

fA  dxfexpCi  r  F(x)]-l}  +  1  .  (C-5) 

Now  performing  the  averages  on  random  variables  ^  and  a^,  we  have,  for  the 
characteristic  function  of  an  individual  component  of  (C-l), 

fl(?)  =  yjjda  p(a)  j'dx^expfi  £ a  F(x)]-^  +  1  ,  (C-6) 

where  p(a)  is  the  probability  density  function  of  random  variable  a^. 
Interchanging  integrals,  (C-6)  becomes 

^(f)  *  Tlldx{fa^F(x)]-1}  +  1  ’  (c~7) 

where  ffl  is  the  characteristic  function  of  amplitude  a^.  Then  from  (C-l), 
since  all  the  individual  random  variables  are  independent,  the  characteristic 
function  of  IK(t)  is 

fi  (fJ  =  [f1(T)3K  •  (C-8 ) 

Finally,  the  characteristic  function  of  total  shot  noise  process  (1)  is,  by 
use  of  discrete  probability  distribution  (A-7)  for  random  variable  K,  given  by 
the  average 
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=  exp[-vT  +  vT  f^F)]  = 

=  exp[vjj"dx{fa[fF(x)]-l}]  .  (C-9 ) 

The  (imprecise)  large  time  interval  T  has  dropped  out  of  the  general 
result  (C-9).  Also,  the  only  parameter  required  about  the  duration  random 
variable  is  its  mean.  The  exact  characteristic  function  f  of 
amplitude  ak  and  the  exact  pulse  shape  F  directly  affect  the  characteristic 
function  of  I(t).  For  =  1  for  all  k,  (C-9)  reduces  to  a  simplified 
version  of  [1,  eq.  1.5-4]. 


Cumulants  of  I(t) 


The  characteristic  function  of  random  amplitude  a^  can  be  expanded  in  a 
power  series 


(C-10) 


where  ua(n)  is  the  n-th  moment  of  a^: 

ua(n)  =  an  =  J  da  an  p(a)  .  (C-ll) 

Then  from  (C-9),  we  develop 
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in  fj(f)  .  vj/dx  {fagF(x)M}  = 

^  •  n  r 

=  »  J  ya(n)  -^ff- Jdx  Fn(x)  , 

n=l 

allowing  for  immediate  identification  of  the  cumulants  of  I(t)  as 


^(n)  =  v  Jua(n)  Jdx  Fn(x)  for  n  >  1;^(0)  =  0  . 
For  J^k  =  1  for  all  k,  this  reduces  to  [1,  eq.  1.5-2]. 

The  normalized  cumulants  of  I(t)  are 


YZ(n) 


£(n)  j  ua(n)  Jdx  Fn(x) 

U^iT  (»I?  '  £a(2)  Jdx  F2(xlT 


These  quantities  tend  to  zero  rapidly  for  >>  1;  see  also  [1,  eq. 
Thus  vj  has  a  pronounced  effect  on  how  Gaussian  I ( t )  is. 


Behavior  of  characteristic  function  f^Cf)  at  f=  *<*> 

If  pulse  F(x)  is  non-zero  only  over(Xp  x^,we  have 

*2 

fdx[fa[fF(x)]-l}  =  \  dx  £faGr F(x)]-l}  . 

J  xx 

Now  if  random  variable  ak  has  a  characteristic  function  fa  with  the 
that 

fa(±*)  =  0 


( C  —12 ) 

(C— 13 ) 

(C-14) 

1.6-3]. 

(C— 15 ) 

property 

(C-16) 
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then 


( C  — 1 5 ) — » 


dx  l  0-l] 


-(x£  -  Xj)  as  %-* 


(C— 17 ) 


in  which  case  (C-9)  yields 

=  exp[-vj(x2  -  xx )]  .  ( C— 18 ) 


If  pulse  extent  x2  -  x^  is  infinite,  as  for  the  Gaussian  or  exponential 


pulses. 


F(x)  =  exp(-x^)  or  exp(-x)U(x) 


( C  —19 ) 


then  (C-18)  is  zero.  On  the  other  hand,  if  x2  -  x^  is  finite,  as  for  circular 
pulse 

r, 

i 


F  ( x )  = 


(l-x2j^for  [x|  < 


0  for  Jxj  >  1 


(C— 20) 


then 


f =  exp[-vj2]  >  0  for  circular  pulse.  ( C  — 2 1 ) 


This  non-zero  characteristic  function  value  corresponds  to  an  impulse  at  the 
origin  of  probability  density  function  Pj,  with  area  (C-21).  Physically, 
this  means  that  there  are  occasionally  regions  of  the  t-scale  where  no^  pulses 
overlap,  and  there  I(t)  =  0.  The  probability  of  this  happening  is,  generally. 


P0  =  Prob  { I (t )  =  oj  =  fj(±<o)  =  exp[-v£(x2  -  x^]  .  (C-22) 
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On  the  other  hand,  for  the  Gaussian  or  exponential  pulses  cited  in 
(C-19),  Xg  -  x^  =  +4  ,  and  fj(±«>)  =  0,  meaning  that  there  is  no  impulse 
at  the  origin  of  probability  density  function  Pj.  Physically,  the  infinite 
tails  (even  if  single-sided,  as  for  the  exponential  pulse)  disallow  I(t)  ever 
from  becoming  zero. 

Cumulants  of  Continuous  Portion  of  Pp 

The  impulse  at  the  origin  means  that  probability  density  function  Pj 
and  cumulative  distribution  function  Pj  can  be  expressed  respectively  as 

P^u)  =  PQ  S(u)  +  pc(u)  , 

u 

Pj (u)  -  PQ  +  S  dt  Pc(t)  for  u  >  0,  (C-23 ) 

where  Pc(u)  is  a  continuous  function  of  u,  with  area  1-PQ.  The 
characteristic  function  relation  corresponding  to  (C-23)  is 

fj(?)  =  P0  +  fc(?)  ,  (C-24) 

and  the  moments  are  related  according  to 

(5l(°)  -  Pq  for  n=0) 

u(n)  J  V  *  (C-25) 

[  Wj(n)  for  n^lj 

The  cumulants  of  f  or  p  can  then  be  found  from  these  moments  (C-25),  by 
recursive  relations;  see  [4]  or  [6].  This  procedure  is  necessary  to  get 
accurate  series  expansions  for  the  probability  density  function  pc  and  its 
cumulative  distribution  function,  without  having  to  approximate  a  delta 
function. 
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In  the  case  where  pulse  extent  x2  -  x^  is  finite,  it  is  possible  to 
find  the  average  number  of  overlapping  pulses  at  any  one  time  instant;  this 
statistic,  denoted  byTp  is  called  the  overlap  factor.  In  order  to 
determine  it  in  a  simple  fashion,  we  concoct  a  very  special  shot  noise 
process:  let 

a^  =  1  for  all  k  , 

F(x)  =1  for  x1  <  x  <  x2  .  (C— 26 ) 

Then  I(t)  is  a  step  function  with  amplitudes  limited  to  the  values 

0,  1,  2,  ...  .  Then  obviously,  the  average  number  of  overlapping  pulses  at 

one  time  instant  is  just 

*i  =  KIT  =  vj  ua(l)  j  dx  F(x)  .  v  J(x2  -  xx)  ,  (C— 27 ) 

upon  use  of  (C-13)  with  n=l  and  (C-26).  If  we  let 

r  =  J(x2  -  Xj)  (C  —28 ) 

denote  the  average  pulse  duration,  we  have  the  overlap  factor  in  the  form 

Kj  =  v  U  .  (C  — 29 ) 

For  the  Gaussian  or  exponential  pulses  in  (C-19),  we  have 
x2  -  x^  =  +«6 ,  giving  L  =  =  +°«.  This  is  in  fact  true,  since 

all  the  infinite  tails  overlap;  however,  it  is  not  then  an  informative 
statistic. 
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Closed  Form  Characteristic  Function  Examples 


There  are  a  couple  of  examples  of  the  circular  pulse  shape  F  and 
amplitude  characteristic  function  fa,  where  (C-9)  can  be  evaluated  in  closed 

a 

form.  This  furnishes  an  alternative  to  the  moment  approach  [4]  used  here. 


Consider  the  circular  pulse  in  (C-20);  then  the  integral  in  (C-9)  is 
(using  (C-l 5 ) ) 

do  coso  fa[?  coso]  -2  ,  (C-30) 

which  holds  for  any  characteristic  function  fa.  Now  first  let  the 

a 

probability  density  function  of  ak  be  exponential: 

P(a)  =  ~exp^-^  U(a)  »  V?)  -  U  -  •  (C-31) 

Substitution  in  (C-30)  yields 


2 


n 


do  coso 

l-ifua  cose 


(C— 32 ) 


But  we  know  that 

it 

2  nr  do  coso 

1-z  coso 


arc  cos(-z)  , 


( C  —33 ) 


via  [5,  eqs.  2.554  2  and  2.553  3].  Then  letting  z  =  i?yfl  and  using 
[2,  eqs.  4.4.2  with  4.4.26],  (C-32)  becomes 
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-2  jji(s-ii  J)  +  i ir (s  — 1 ) 


wa?s 


-2,  with  s 


-  (>  *  4ft- 


( C  —34 ) 


Combining  these  results  in  (C-9),  the  closed  form  characteristic  function  is 

fj(?)  =  exP  j~  +^n(s-pa|)  -  i  j(s-l)]J  ,  (C— 35 ) 

which  holds  for  a  circular  pulse  F  and  an  exponential  probability  density 
function  p(a). 


The  second  example  is  the  one  considered  in  detail  here,  namely  the 
Rayleigh  probability  density  function  p(a)  given  in  (3).  First  substituting 
(C-30)  in  (C-9),  we  have  characteristic  function 


fj(f)  =  exp[2v4(J(?)-l)] 


where  integral  J(^)  is  defined  as 


J(f)  = 


de  cose  f  JJ  cose] 


(C— 36 ) 


(C-37) 


For  Rayleigh  probability  density  function  (3),  (C-37)  can  be  expressed  as 
follows: 


J(f) 


IT 

\  de  cose 
0 


r 


da  exp(iaj  cose) 


(C-38) 


Transform  to  rectangular  coordinates  according  to  a  cos  e  =  aa  x, 
a  sin  e  =  aa  y,  and  obtain 


J(2) 


dx  dy 


(x 


V)5 


6 


exp(ifaax 


(C-39) 
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But  the  integral  on  y  here  is,  via  y  =  x  u  ,  equal  to 

?  r  ITi:  exp(-  f  A)  -  ?  exp(r)  Ko(r)  •  <CJ>°) 

0  u'(l+u) 

the  latter  by  means  of  [5,  eq.  3.364  3].  Thus  (C-39)  becomes 
J(£)  =  J  dxx  exp^ifaax-  ^  exp (^j  KQ 


At  this  point,  we  have  two  alternatives.  First,  (C-41)  could  be 
efficiently  evaluated  for  all  S’  via  an  FFT;  the  decay  of  the  integrand  is 
according  to  exp(-u  )  for  large  u.  Secondly,  J(/)  can  be  expressed  in  a 
closed  form  in  terms  of  a  hypergeometric  function;  specifically 

J(f)  =  2F2(1,  1;  £;  -2b2)  + 

+  i(|J^  b  exp(-b2)  [IQ(b2)  -  I^b2)]  ,  (C -42) 

where  b  =  cra£/2.  The  upper  line  follows  from  [5,  eq.  6.755  6],  while  the 
lower  line  used  [5,  eq.  6.755  9]  with  an  application  of  partial  derivative 
3/sa  to  both  sides.  The  characteristic  function  fj  is  finally  obtained  by 
employing  (C-42)  in  (C-36). 

Still  another  alternative  is  afforded  by  use  of  the  closed  form  for  the 
characteristic  function  of  the  Rayleigh  probability  density  function,  as  given 
in  [7,  eq.  6]. 


du  exp(i  2  of  u)  u  exp 


4 


t) 


(C-41) 
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Moments  for  Some  Particular  Pulse  Shapes  F 


The  moments  of  pulse  shape  F  were  encountered  in  evaluating  the 
cumulants  %j(n)  of  shot  noise  I(t),  according  to  (24)  or  (C-13)  as 


jp(n)  =  Jdx  Fn(x)  for  n  >  0 


For  circular  pulse  (C-20),  [5,  eq.  3.621  1]  yields  moments 


p  11 

up(n)  =  T  dx  (l-x2)*  =  2  d©  (cos©) 


-1 


a  result  already  quoted  in  (29). 


?n+lp2/n+1N 

n+1  2  r  It1) 

-  p(n+2)  ’ 


More  generally,  for 


F(x)  = 


(l-x2)a  for  |x|  <  1 


0  for  |x|  >  1 
[5,  eq.  3.621  1]  yields,  with  a  trigonometric  substitution. 


/„\  „2na+l  r(na+l) 
uF(n)  =  2  r(W 2)  • 


For 


F(x)  = 


}(cos  x)a  for  |xj  <  j 
0  for  lx)  >  | 


[5,  eq.  3.621  1]  yields  directly 

/na+1' 


P2(- 

-  2naTrfei7 


(C-43) 


(C-44 ) 


(C -45 ) 


(C-46) 


(C-47 ) 


( C  —48 ) 
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For 

F(x)  =  xa  exp(-x)  U(x)  , 

»F<n>  fc-49* 

n 

while  for 

F (x)  =  xa  exp(-x2/2)  U(x)  , 


Both  relations  follow  directly  from  the  definition  of  the  P  function. 

Some  Probability  Density  Functions  for  Amplitude  ak 
For  probability  density  function 

PU)  ■  ula>  ■  (c-51) 

we  have  characteristic  function 

fa(f)  =  (1-ifa)^'1  CC— 52) 

with  moments 

Ua(n)  =  (Y+l)n  <*"  for  n  >  0  (C-53) 

and  cumulants 

^(n)  =  (n-l)l  (y+1  )  «n*  for  n  )  1  .  (C-54) 

This  example  subsumes  the  exponential  probability  density  function,  upon 
setting  y  =  0. 
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For  probability  density  function 

nm  2(a/g)r  exp(-a2/g2) 


0  .  (C-56) 

This  example  subsumes  respectively  the  one-sided  Gaussian  for  y=°>  the 
Rayleigh  for  y=1»  and  the  Maxwell  probability  density  functions  for  y=2.  The 
result  in  (29)  follows  immediately  by  setting  y=1»  a =  2'* a  . 

a 

Convergence  of  Series  for  j?n  f|(f) 

A  power  series  expansion  for  j^n  fj(f)  was  developed  in  (C-12),  namely. 

In  f  j(f)  =  ua(n)  uF(n)  ;  (C-57) 

n=l 

here  we  employed  (C-43).  Since  the  moments  in  (C-57)  can  be  easily  evaluated 
via  recursion,  according  to  results  in  the  above  two  subsections,  it  might  be 
thought  that  (C-57)  could  be  employed  to  evaluate  the  characteristic  function 
of  I(t)  directly,  without  recourse  to  the  more  difficult  approaches  required 
in  (C-9 )  or  (C-35)  or  (C-36)-(C-42) . 

To  see  the  drawbacks  of  this  approach,  consider  first  a  circular  pulse  F 
and  a  generalized  exponential  probability  density  function  p(a)  as  in  (C-51); 
then  a  combination  of  (C-44)  and  (C-53)  yields,  for  the  n-th  term  of  the  sum 
in  (C-57), 


we  have  moments 


Ua(H) 


rf^  „ 


for  n  > 
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.2<rl)/gi) 

“  ni(n+l): 


(ifa2)n. 


(C-58) 


Then  ratio 

n — =  _  ip+y  )(n+y-l)  y2a2^  _  j2a2  as  n-»  +•«,  (C-59) 

n-2  n  -1 


regardless  of  the  value  of  y.  Therefore 


n-1 


^  Jf  j a  as  n-»  +  *«  , 


(C-60) 


meaning  that  series  (C— 57 )  only  converges  for  \f\<  1/a.  So  (C-57)  is  not  a 
viable  approach  for  the  calculation  of  the  characteristic  function  in  this 
case. 


As  a  second  example,  we  consider  the  circular  pulse  F  with  the 
generalized  Rayleigh  probability  density  function  in  (C-55).  Combination  of 
( C  —44 )  with  ( C— 56 )  yields  for  the  n-th  term  of  series  (C-57), 


T 


n 


n!  (n+l)  1 


(i?a2)n  . 


(C  —61 ) 


Then  ratio 


_ ?_  =  .  H.tr,-.1  „  _  I2?2 

Tn-2  n2-l  2  2" 


as  n-»  +  , 


regardless  of  y.  Therefore 


I 


as  , 


(C— 62 ) 


(C-63 ) 
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meaning  that  series  (C-57)  converges  for  all/*.  However,  direct  numerical 
evaluation  of  (C-57)  via  (C-61)  and  (C-62)  loses  all  its  significant  digits 
for  large  long  before  £n  fj(J)  reaches  its  final  value  of  -2J?v  +  iO,  due 
to  the  alternating  character  of  the  series.  So  (C-57)  is  not  a  useful 
approach  for  evaluation  of  the  characteristic  function,  except  for  small  f. 
By  contrast,  the  series  expansion  technique  employed  in  [4]  uses  the  moments 
to  directly  estimate  the  desired  probability  density  function  and  cumulative 
distribution  function  of  interest,  for  large  arguments  as  well  as  small. 


51/52 
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APPENDIX  D.  PROGRAMS  FOR  CUMULATIVE  DISTRIBUTION  FUNCTION 
AND. PROBABILITY  DENSITY  FUNCTION 

The  programs  used  here  to  evaluate  the  cumulative  distribution  function 
and  probability  density  function  of  shot  noise  are  listed  below.  The  n-th 
coefficient  in  a  generalized  Laguerre  expansion  of  orthonormal  polynomials  is 
denoted  by  bn  and  is  plotted  in  figure  D-l  for  n  =  0(1)70.  It  is  seen  to 
oscillate  and  decay  with  n  until  n  =  32,  at  which  point  round-off  error 
becomes  important;  however,  by  this  time,  |b  |  has  decayed  below  the  IE-5 
level.  The  round-off  error  is  so  dominant  beyond  n  =  35,  that  no  useful 
results  for  bn  can  be  obtained  then.  The  particular  parameter  values  (a,  b) 
used  for  the  Laguerre  weighting  are  indicated  in  the  listings. 
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10 
20 
30 
40 
50 
60 
70 
30 
90 
100 
110 
120 
130 
140 
150 
166 
170 
130 
190 
200 
210 
220 
230 
240 
250 
260 
270 
280 
290 
300 
310 
320 
330 
340 
350 
360 
370 
330 
390 
400 
410 
420 
430 
440 
450 
460 
470 
480 
490 
500 
510 
520 
530 
540 
550 
560 
570 
'  580 
590 
600 
610 
620 


STEP  PLUS  CONTINUOUS  PART  OF  SHOT  NOISE  CDF,  Pc<u>, 

VIA  GENERALIZED  LAGUERRE  EXPANSION  AND  MOMENTS 

M=70  !  MAXIMUM  ORDER  OF  APPROXIMATION;  NUMBER  OF  MOMENTS  REQUIRED 


DOUBLE  M, I , N, K  ! 

REDIM  Mom(0:M>,A<0:M>,L<0sM> 

REAL  Mom < 0 : 1 00 > , A < 0 : 1 00 > , L < 0 : 1 00 > , 
CALL  Moment s ( M, P@, Mom <*) >  ! 

Center=Mom<l)^Mom<0>  ! 

R2=Mom<2>''Mom(0>-Center*Center  ! 
Rms=SQR  <  R2  >  ! 

A 1  pha0=Cent  er*Cent er^R2- 1 .  ! 

Beta0=R2^Cent er  ! 


INTEGERS  <  2 A 3 1  =  2,147,483,648 
Y< 1 : 21 > 

P 0  IS  STEP  AT  ORIGIN 
CENTER  OF  pc<u> 

MEAN  SQUARE  SPREAD  OF  pc(u> 

RMS  SPREAD  OF  pc<u> 

THE  CHOICES  A 1  pha=A1  pha@  AND 
Beta=Beta@  WOULD  MAKE  A<l)=fl<2> 


=  0 


A 1 pha= . 74 


Beta=2. 1 


CALL  Coeff 1 d_y i a_mom<M, A1 pha. Beta, MomC*) , A<*> >  !  DIRECT  MOMENT 

CALL  Coeff 1 r~gi a_mom<M, A1 pha, Beta, Mom<*) , A<*> )  !  RECURSIVE  MOMENT 

PRINT  "Center  =  ";Center 
PRINT  "Rms  =  "  ;  R  m  s 
Al=Alpha+l. 

01=1. /fll 

F 1  =  1 .  /'FNGamma<  A 1 ) 

DATA  . 001 , . 002, . @05, . 01 , . 02, . 05, . 1 , . 2, . 3, . 4, . 5 
DATA  . 6, . 7, . 3, . 9, . 95, . 98, . 99, . 995, ’ 998, . 999' 

READ  Y<*> 

FOR  1=1  TO  21 
Y<I>=FNInvphi < Y  <  I  >  > 

NEXT  I 
Y 1  =Y  <  1 ) 

Y2  =  Y  <  21) 

INPUT  "ORDER  AND  L IM I TS : " , N , U 1 , U2 
PRINT  "ORDER  AND  L I M I TS : " , N ; U 1 ; U2 
Du=<U2-Ul>/'100. 

PLOTTER  IS  "GRAPHICS" 

GRAPHICS  ON 

WINDOW  U 1 , U2, Y 1 , Y2 

FOR  U=U1  TO  U2  STEP  <U2-Ul)*.l 

MOVE  U , Y 1 

DRAW  U,Y2 

NEXT  U 

FOR  1=1  TO  21 
MOVE  U 1 , Y < I ) 

DRAW  U 2 , Y ( I ) 

NEXT  I 
PENUP 

FOR  1=1  TO  100 
U=Ul+Du*I 
T  =  U  /  B  e  t  a 

CALL  Laguerre<N-l , A 1 , T , L< *> > 

Sum=A<0)*FNFll<Al,T)*Ol 
FOR  K= 1  TO  N 


Sum  =  Sum  +  fl<K)*L<K-l)<''K 
NEXT  K 

P=P0+Fl*EXP<-T+Al*LOG<T)>#Sum  '  !  PROBABILITY  THAT  RV  <  U 
IF  P>0.  AND  P<1.  THEN  570 
PENUP 
GOTO  580 

PLOT  U, FNInvphi  <.  P ) 

NEXT  I 
PENUP 
GOTO  300 

END 

I 
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639  B  E  F  FNInMphi(X)  !  INVPHIiX)  Mia  26.2.23  with  mo d i f i c  at  i o  n 

648  D=X-. 5 

650  IF  HBS ( D ) > . 0 1  THEN  680 

660  P=2. 50662827463*0* < 1 . +D*D*1 . 04719755120) 

670  RETURN  P 

680  P=X 

690  IF  X > . 5  THEN  P=l.-X 

700  P=SQR<-2. *LQG  <  P ) ) 

710  T= 1 . +P*< 1. 432783+P*< . 1 39269+P* . 00 1 308 > ) 

720  P  =  P- <  2 . 515517+P*< . 802S53+P*. 010328) )  /  T 
730  IF  X< . 5  THEN  P=-P 

740  RETURN  P 

750  FNEND 

769  ! 

770  DEF  FNGammaCX)  !  GammaCX)  via  HART,  page-  282,  #5243 

780  DOUBLE  N ,  K 

790  N  =  I  N  T ( X  ) 

800  R=X-N 

810  IF  N>0  OR  R<  >0 .  THEN  840 

820  PRINT  "FNGamma(X)  IS  NOT  DEFINED  FOR  X  =  "  ;  X 

830  STOP 

840  IF  R >0 .  THEN  870 

850  Gamma2=l. 

860  GOTO  940 

870  P=439. 33@444060025676+R*<50. 1 086937529709530+R*6 . 74495072459252899 ) 

880  P=8 762. 7 10297352 14896+R*< 2003. 52740 1307279 12+R*P) 

890  P=42353. 6895097440396+R*<20886. 36 1 7892698874+R*P ) 

900  Q=4 99. 02852662 143904S-R*( 189. 498234 1 570280 1 6-R* < 23 . 08 1 55 1 524580 1 25-R ) ) 

910  Q=9940. 307415082770 90-R*< 1528. 607273 7 79 5220 2 +R*Q) 

920  Q=42353. 6895097440900+R*< 2980. 38533 0925 66499- R-*-Q) 

930  Gamma2=P^Q  !  Gamma<2+R>  for  0  <  R  <  1 

940  IF  N >2  THEN  980 

950  IF  N<2  THEN  1030 

960  Gamma=Gamma2 

970  RETURN  Gamma 

980  Gamma=Gamma2 

990  FOR  K= 1  TO  N-2 

1000  Gamma=Gamma*(X-K) 

1010  NEXT  K 

1020  RETURN  Gamma 

1030  R= 1 . 

1040  FOR  K=0  TO  1-N 

1050  R=R*<X+K) 

1060  NEXT  K 

1070  Gamma=Gamma2/R 

1080  RETURN  Gamma 

1090  FNEND 

1100  ! 
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1110  DEF  FNF 1 1 <  ft  1 , X)  !  1F1 < 1 ; fil+1 ; X> 

1120  DOUELE  K 

1130  T=S=1 . 

1140  FOR  K= 1  TO  20© 

1150  T=T*X/<fll+K> 

1160  S  =  S  +  T 

1170  IF  T<=1.E-17*S  THEN  RETURN  S 
1180  NEXT  K 

1190  PRINT  "200  TERMS  IN  FNF 1 1  RT“ ; fll ; X 
1200  RETURN  S 

1210  FNENB 

1220  ! 

1230  SUE  Laguerre < DOUELE  H,RERL  fll pha, X,L(*) )  !  Lri\al  pha(X> 

1240  DOUBLE  K 

1250  Rl=fl1 pha-1 . 

1260  L  <  0 )  =  1 . 

1270  L<  1 ) =A1 pha+ 1 . -X 

1280  FOR  K=2  TO  N 

1290  L<K)  =  <<K+K+FI1-X>*L<K-1  >-<K+fln*L<K-2>VK 

1300  NEXT  K 

1310  SUBEND 

1320  ! 

1330  SUB  Momnt  via  cumnt< DOUBLE  M,REflL  Cum < , Mom ( * > ) 

1340  DOUBLE  K,J 

1350  REAL  Mom0 

1360  Mom<0>=Mom0  =  EXP<Cum<0>  > 

1370  FOR  K= 1  TO  M 

1380  T=1 . 

1390  S=Cum<K)*Mom0 

1400  FOR  J=1  TO  K-l 

1410  T=T*<K-JW 

1420  S=S+T*Cum<K-J>*Mom(J) 

1430  NEXT  J 

1440  MomOO=S 

1450  NEXT  K 

1460  SUBEND 

1470  ! 
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1480  SUE  Coef f 1 d_vi a_Bom( DOUBLE  M , REAL  A1 pha. Beta, Mom(*>  ,  A<*> > 

1490  ALLOCATE  B<0:M> 

1500  DOUBLE  K ,  K 1 ,  J ,  M  x 

1510  T=l. 

1520  FOR  K= 1  TO  M 

1530  T=T*< H 1 pha+K ) *Eet a 

1540  MomOO=MomaOxT  !  NORMALIZED  MOMENTS,  RELATIVE  TO  Alpha  AND  Bet 
1550  NEXT  K 

1560  Q= 1 . 

1570  A<0)=B<0>=Mom<0) 

1580  FOR  K= 1  TO  M 

1590  K 1 =K+1 

1600  T= 1 . 

1610  S=Mom<0) 

1620  FOR  J=1  TO  K 

1630  T=T*< J-Kl )xj 

1640  S=S  +  T*Morn  <  J  !> 

1650  NEXT  J 

1660  Q=Q*<Alpha+K)xK 

1670  A  <  K ) =S 

1680  B<IO=S*SQR<Q) 

1690  NEXT  K 

1700  Mx=Mx+10 

1710  IF  Mx<  M  THEN  1700 

1720  Threshol d=-7. 

1730  T2=Thr esho 1 d*2 . 

1740  V=10. "Threshol d 

1750  GINIT 

1760  PLOTTER  IS  "GRAPHICS" 

1770  GRAPHICS  ON 

1780  WINDOW  0. ,FLT<Mx>,T2,0. 

1790  LINE  TYPE  3 

1800  FOR  J=0  TO  Mx  STEP  10 

1810  MOVE  J , T2 

1820  DRAW  J,0. 

1830  NEXT  J 

1340  FOR  J=T2  TO  0 

1850  MOVE  0. , J 

I860  DRAW  Mx , J 

1870  NEXT  J 

1880  PENUP 

1890  LINE  TYPE  1 

1900  IMAGE  4D , 2 ( 4X , M . 1 7DE > 

1910  PRINT  "  K  BOO  Sum" 

1920  Sum=0. 

1930  FOR  K=0  TO  M 

1940  B=B<K) 

1950  Sum=Sum+B*B 

I960  PRINT  USING  1 900; K , B , Sum 

1970  IF  B< V  THEN  2000 

1980  Y  =  LGT  (B) 

1990  GOTO  2040 

2000  IF  B >-V  THEN  2030 

2010  Y=T2-LGT  < -B ) 

2020  GOTO  2040 

2030  Y=Threshold 

2040  PLOT  K , Y 

2050  NEXT  K 

2060  PENUP 

2070  SUBEND 

2080  ! 
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2090  SUB  Coef  f  1  r_u  i  ajioni  <  DOUBLE  M ,  REAL  A  1  p  h  a ,  Beta,  Morn  (  *  }  ,  A  C  *  > ) 

2100  ALLOCATE  B<0:M> 

2110  DOUBLE  K , K 1 , J , M x 

2120  T=l. 

2130  FOR  K= 1  TO  M 

2140  T  =  T*<Alpha+IO*Beta 

2150  MomOO=Mom<K>/T  !  NORMALIZED  MOMENTS,  RELATIVE  TO  Alpha  AND  Bet 
2160  NEXT  K 

2170  Q=  1  . 

2180  A0  =  A <  0 )  =B  <  0  >  =  M o rn  <  0  > 

2190  FOR  K= 1  TO  M 

2200  K 1 =  K+ 1 

2210  T=  1 . 

2220  S=Mom<K)-A0 

2230  FOR  J=1  TO  K-l 

2240  T=T*(J-K1 )/J 

2250  S=S-T*A< J> 

2260  NEXT  J 

2270  IF  K  MOD  2=1  THEN  S=-S 

2280  Q=Q*<A1 pha+K>/K 

2290  A ( K ) =S 

2300  B<K>=S*SQR<Q> 

2310  NEXT  K 

2320  Mx=Mx+10 

2330  IF  M x < M  THEN  2320 

2340  Thresho 1 d=-7 . 

2350  T2=Threshold*2. 

2360  V= 1 0 . AThresho 1 d 

2370  GINIT 

2380  PLOTTER  IS  "GRAPHICS" 

2390  GRAPHICS  ON 

2400  WINDOW  0. ,FLT<Mx),T2,0. 

2410  LINE  TYPE  3 

2420  FOR  J=0  TO  Mx  STEP  10 

2430  MOVE  J , T2 

2440  DRAW  J,0. 

2450  NEXT  J 

2460  FOR  J=T2  TO  0 

2470  MOVE  0.,J 

2480  DRAW  Mx , J 

2490  NEXT  J 

2500  PENUP 

2510  LINE  TYPE  1 

2520  IMAGE  4D , 2 < 4X , M . 1 7DE > 

2530  PRINT  "  K  BOO  Sum" 

2540  Sum=0. 

2550  FOR  K=0  TO  M 

2560  B=B(K> 

2570  Sum=Sum+B*E 

2580  PRINT  USING  2520; K,B, Sum 

2590  IF  B< V  THEN  2620 

2600  Y  =  LGT  <B) 

2610  GOTO  2660 

2620  IF  E >-V  THEN  2650 

2630  Y  =  T2-LGT  <-E) 

2640  GOTO  2660 

2650  Y=Threshold 

2660  PLOT  K, Y 

2670  NEXT  K 

2680  PENUP 

2690  SUBEND  59  N 

2700  ! 
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2710 

2720 

2730 

2740 

2750 

2760 

2770 

2780 

2790 

2800 

2810 

2820 

2830 

2840 

2850 

2360 


SUE  M oments< DOUBLE 
Over  1 ap=6 . 2 
S i gmaa= 1 . 
F'0=EXP<-Over1  ap) 
ALLOCATE  Cum<0:M> 


M , REAL  P@, Mom<*) )  !  SHOT  NOISE 

!  AV.  NO.  PULSES/'SEC  *  AVERAGE  PULSE  DURATION 
!  PARAMETER  OF  RAYLEIGH  AMPLITUDE  PDF 
!  PROBABILITY  OF  ZERO  AMPLITUDE  OF  SHOT  NOISE 
!  ARRAY  FOR  CUMULANTS 


DOUBLE  K 
S  =  S  i  gmaa*S i  grnaa 
Cum  <  0 ) =0 . 

Cum<  1  )=Over1  ap*Si  gmaa*.  25*PI*SQR< .  5*PI  > 
Cum<2)=0ver1ap*S*4.-"3. 

FOR  K=3  TO  M 

Cum<K)=Cum<K-2>*S*K*K/<K+l) 


NEXT  K 

CALL  Momnt_v i a_c urnnt  <  M , Cum <  *  > , Mom  <  *  )  ) 
Mom<0)=Mom<0)-P0 


SUBEND 


10 
20 
30 
40 
50 
60 
70 
80 
90 
100 
1  10 
120 
130 
140 
150 
160 
170 
180 
190 
200 
210 
220 
230 
240 
250 
260 
270 
280 
290 
300 
310 
320 
330 
340 
350 
360 
370 
380 
390 
400 
410 
420 
430 


CONTINUOUS  PART  OF  SHOT  NOISE  PDF,  pc<u>,  VIA 
GENERALIZED  LAGUERRE  EXPANSION  AND  MOMENTS 

M=90  !  MAXIMUM  ORDER  OF  APPROXIMATION;  NUMBER  OF  MOMENTS  REQUIRED 

DOUBLE  M , I ,  N ,  K  !  INTEGERS  <  2A31  =  2,147,483,648 


R E D I M  Mom<0:M),A<0:M),L<0sM) 

REAL  Mom < @ : 100) , A<0; 100) , L<0: 100) 
CALL  Moment s <M, P0, Mom <*) )  ! 

Center=Mom(  1  )/Morn<0)  ! 

R2=Mom<2)/Mom<8) -Cent er*Cen ter  ! 
Rm*-SQR<R2)  ! 

A 1 pha8  =  Cent  er*Cent  er/R2-l .  ! 

BetaO=R2^Center  ! 

A 1 pha= . 74 
Bet  a=2 . 1 
CALL 
CALL 
PRINT 

PRINT  "  Rms  =";Rms 
F 1  =  1 . ✓<  Beta*FNGamma< A1  pha+1 .  >) 

INPUT  "ORDER  AND  LIMITS; ", N, U1 , U2 
PRINT  "ORDER  AND  LI M I TS : " , N ; U 1  $ U2 
Du=<U2-U 1)^100. 

H  =  4 .  /  (.  U2-U 1 ) 

PLOTTER  IS  "GRAPHICS" 

GRAPHICS  ON 
WINDOW  U1,U2,-H*. 1,H 
GRID  <U2-U1 )*. 1 , H*.  1 
PLOT  0. , 0. 

FOR  1=1  TO  100 
U=Ul+Du*I 


P 0  IS  STEP  AT  ORIGIN 
CENTER  OF  pc<u> 

MEAN  SQUARE  SPREAD  OF  pcCu) 

RMS  SPREAD  OF  pc<u) 

THE  CHOICES  Alpha® AlphaO  AND 
Bet a=Bet aO  WOULD  MAKE  A<1)=A<2)=0 


Coet'f  1  d_vi  a_mom<M,  A1  pha,  Beta,  Mom<*)  ,  AO)  )  !  DIRECT  MOMENTS 

Coeff 1 r_vi a~mom<M, A1 pha, Beta, Mom<*> , fl<*) )  !  RECURSIVE  MOMENTS 

"Center  =  Center 


T=U^Be ' a 

CALL  Laguerre<N, A1 pha, T, L<») ) 

Sum=A<0) 

FOR  K= 1  TO  N 
Sum=Sum+ A  <K)*L<K) 

NEXT  K 

P=Fl*EXP(-T+Rlpha*LQG<T))*Sum  !  PDF  OF  RV  AT  U 

PLOT  U,P 

NEXT  I 

PENUP 

GOTO  200 

END 


i 
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